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ABSTRACT 



We zoom in on the internal structure of the low-frequency quasi-periodic oscillation (LF QPO) often observed in black hole binary 
systems to investigate the physical nature of the lack of coherence in this feature. We show the limitations of standard Fourier 
power spectral analysis for following the evolution of the QPO with time and instead use wavelet analysis and a new time-frequency 
technique - Matching Pursuit algorithm - to maximise the resolution with which we can follow the QPO behaviour We use the LF 
QPO seen in a very high state of XTE J1550-564 to illustrate these techniques and show that the best description of the QPO is that it 
is composed of multiple independent oscillations with a distribution of lifetimes but with constant frequency over this duration. This 
rules out models where there is continual frequency modulation, such as multiple blobs spiralling inwards. Instead it favours models 
where the QPO is excited by random turbulence in the flow. 

Key words. Methods: data analysis ~ Techniques: miscellaneous - X-ray: binaries - X-ray: individuals: XTE J1550-564 



1. Introduction 

One of the most outstanding features of the stellar mass black- 
hole binaries (BHB) is their rapid X-ray time variability (e.g. van 
der Klis 1989). This variability consists of aperiodic continuum 
noise spanning a broad range in frequencies, with quasi-periodic 
oscillations (QPOs) superimposed. The 'low frequency QPO' 
(hereafter LFQPO) is the one most often seen, and this is at a 
characteristic frequency which moves from 0.1-10 Hz in a way 
which is correlated with the energy spectrum of the BHB (see 
e.g. the reviews by van der Klis 2006; McClintock & Remillard 
2006; Done et al. 2007). 

The true nature of this LFQPO (or any of the other QPO's) 
is still not well understood, and there are multiple suggestions 
in the literature. The most fundamental frequency is Keplarian, 
but this is very high, at ~ 200 Hz at the last stable orbit for 
a IOMq Schwarzchild black hole (Sunayev 1973). Instead, the 
QPO could be linked to the buildup and decay of 'shots' on the 
longer timescales for magnetic reconnection events (Miyamoto 
& Kitamoto 1989; Negoro et al. 1994), though to form a QPO 
probably requires correlation between the shots via an energy 
reservoir (Vikhlinin et al. 1994). Alternatively, the frequency 
could be set by modes excited in the accretion disc (Nowak 
& Wagoner 1995; Nowak et al. 1997), or by one-armed spi- 
ral waves in the disc (Kato 1989) or Lense-Thirring preces- 
sion (Ipser 1996; Stella & Viertri 1998: Psaltis & Norman 2000; 
Ingram et al. 2009). 

Clearly, observational constraints are required in order to 
distinguish between some of these very dilferent physical mod- 
els. Firstly, the energy dependence of the QPO clearly shows 
that the oscillation does not involve the optically thick accre- 
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tion disc emission but is instead a feature of the higher energy 
Comptonised tail (e.g. Zycki et al. 2007). Therefore all mod- 
els using a modulation of the thin disc are ruled out. Instead, the 
QPO must be connected to the hard X-ray emitting coronal mate- 
rial, though the 'shot noise' models are ruled out by the observed 
r.m.s. -fiux correlation as this cannot be formed by independent 
events (Uttley et al. 2005). 

More constraints come from the detailed shape of the QPO 
signal. The most popular method for studying X-ray variabil- 
ity is the Fourier power spectrum density (PSD; see Vaughan 
et al. 2003 and references therein). This is the square of the 
amplitude of variability at a given frequency, as a function of 
frequency, namely P(f). Power spectra from BHB can be very 
roughly described as band limited noise, with a 'flat top' in vP{v) 
(equal variability power per decade in frequency) i.e. P{v) cc v"' . 
This extends between a low and high frequency break, Vb, below 
which the PDS is P(v) oc v", and v/,, above which the spectrum 
steepens to P{v) oc v"^. However, this broadband noise is better 
described by a sum of 4-5 Lorentzian components (Belloni & 
Hasinger 1990; Psaltis et al. 1999; Nowak 2000; Belloni, Psaltis 
& van der Klis 2002). This has the advantage that LF QPO is 
also well modelled by a Lorentzian, so both QPO and under- 
lying noise can be fitted with the same functions. More impor- 
tantly, Lorentzians also suggest a physical interpretation of the 
power spectral components as these are the natural outcome of a 
damped, driven harmonic oscillator. Each Lorentzian-like struc- 
ture in a PSD may then be related to a process in the accretion 
flow, hence X-ray emission x(t) oc sin(27r/()f) exp(-f/T) where r 
denotes a damping time-scale (Misra & Zdziarski 2008; Ingram 
et al. 2009). 

A Lorentzian has P(f) = AA//[(/ - fof + (Af/2)% de- 
fined by a centroid frequency, /o, a full width at half maximum 
(FWHM) of A/ and amplitude. A ratio of the QPO frequency to 
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Lorentzian FWHM is known as a quality factor, Q - /o / A/. The 
distinguishing feature of the QPO as opposed to the noise com- 
ponents is that its quahty factor can be large. All these paramters 
are correlated in the LF QPO, with the amplitude and quality fac- 
tor increasing along with the increase in centroid frequency (e.g. 
Nowak et al. 1999; Nowak 2000; Pottschmidt et al. 2003). 

It is obviously important to know what sets the (lack of) co- 
herence of the QPO. Are QPOs composed of a set of longer- 
lasting continuous modulations or are they rather fragmented in 
time with frequencies concentrated around the peak QPO fre- 
quency? Is there any correlation between observed oscillations 
or they are preferably excited in a random manner? To answer 
these questions we need to resolve and understand the internal 
structure of detected QPOs i.e. to study the behaviour of the 
QPO on timescales which are short compared to its observed 
broadening. A continuous modulation would then show up as 
a continuous drift in QPO frequency with time, while a series 
of short timescale oscillations will show up as disjoint sections 
where the QPO is on and off. 

One way to do this is using dynamical power spectra (spec- 
trograms; hereafter also referred to as STFT for simplicity). 
While the actual techniques for this can be quite sophisticated, 
(e.g. Cohen 1995; Flandrin 1999), in essence, an observation of 
length T sampled every At (giving a single power spectrum span- 
ning frequencies l/T to 1/(2 A?) in steps of l/T) is spilt into 
segments of length T/N. This gives A^ independent power spec- 
tra spanning frequencies N/T to l/(2Af) but crucially, the reso- 
lution is now lower, at N/T. While this is useful in tracing the 
evolution of the QPO (e.g. Wilms et al. 2001; BaiTet et al. 2005), 
it introduces 'instrumental' frequency broadening from the win- 
dowing of the data which prevents us following the detailed be- 
haviour of the QPO on the required timescales. 

The real problem with such Fourier analysis techniques is 
that they decompose the lightcurve onto a basis set of sinusoid 
functions. These have frequency /, with resolution A/ = N/T, 
but exist everywhere in time across the duration of the observa- 
tion tdur - T/N. In the same way that the Heisenburg uncertainty 
principle sets a limit to the measurement of momentum and lo- 
cation of a particle, namely Ap^ Ax > h/ln where h is a Planck 
constant, there is a Heisenberg-Gabor uncertainty principle set- 
ting the limiting frequency resolution for time-series analysis. 
This states that we cannot determine both the frequency and time 
location of a power spectral feature with infinite accuracy (e.g. 
Flandrin 1999). 

Instead, if the QPO is really a short-lived signal, we will gain 
in resolution and get closer to the theoretical limit by using a set 
of basis functions which match the underlying physical shape of 
the QPO. This is the idea behind wavelet analysis. For example, 
one particular basis function shape is the Morlet wavelet, which 
is a sinusoid of frequency /, modulated in amplitude by a gaus- 
sian envelope such that it lasts only for a duration tjur. The prod- 
uct /xrj„r is set at a constant, fixing the number of cycles seen in 
the basis function. The lightcurve is then decomposed on these 
basis functions, calculated over a set of frequencies so that the 
basis function shape is maintained (i.e. that the oscillation con- 
sists of 4 cycles, so low frequencies have longer t^ur than high 
frequencies). The resolution adjusts with the frequency, making 
this a more sensitive technique to follow short duration signals 
(e.g. Farge 1992; Addison 2005). 

Quite quickly, wavelet transforms gained huge popularity 
within the astronomical community (e.g. Szatmary, Vinko & Gal 
1994; Frick et al. 1997; Aschwanden et al. 1998; BaiTeiro & 
Hobson 2001 ; Freeman et al. 2002; Irastorza et al. 2003) whereas 
its application in X-ray astronomy was not so rapid. Scargle et al. 



(1993) used wavelets to examine QPO and very low-frequency 
noise in Sco X-1. Steiman-Cameron et al. (1997) supplemented 
Fourier detection of quasi-periodic oscillations in optical light 
curves of GX 339^ by wavelets while Liszka et al. (2000) anal- 
ysed the ROSAT light curve of NGC 5548 in short time-scales. 
Recent studies of QPOs in X-ray black-holes systems can be 
found e.g. in Lachowicz & Czerny (2005), Espaillat et al. (2008) 
and Gupta et al. (2009). 

However, the problem is that the basis functions chosen in 
wavelet analysis may not be appropriate. For example, we as- 
sumed above that these functions have a fixed shape, lasting for 
a fixed number of oscillations at all frequencies. In practice, this 
may not be the best description of the QPO. Perhaps the QPO is 
made from a set of signals which have a distribution of durations, 
where / x frf„,. is not constant. We will maximise the resolution 
with which we can look at the QPO if and only if we use basis 
functions which best match its shape. 

To do this we propose the application of the Matching 
Pursuit algorithm (MP). This is an iterative method for signal de- 
composition which aims at retrieving the maximum possible the- 
oretical resolution by deriving the basis functions from the sig- 
nal itself. We specifically use this as an extension of the wavelet 
technique by setting the MP basis functions as Gaussian am- 
plitude modulated sinusoids as before, but allowing the product 
/ X tijur to be a free parameter (Gabor atoms). 

We apply both wavelet and Matching Pursuit analysis to 
zoom in on the detailed structure of the LF QPO. A technical 
description of these time-frequency techniques we provide in 
the Appendix. We describe the data used in Section |2l This is 
a RXTE/PCA observation of the BHB XTE J 1550-564 which 
which displays a strong and narrow QPO at ~4 Hz. Section 
[3] shows how wavelets and MP give more information than 
the standard dynamical power spectra. We study wavelet and 
Matching Pursuit feedback to different QPO models in Section|4] 
whereas in Section |5] we discuss the new physical insights this 
gives about the QPO mechanism. 

2. Source selection and data reduction 

XTE J 1550-564 is X-ray nova and black-hole binary discovered 
on September 6, 1998 (Smith et al. 1998) by All Sky Monitor 
on-board the Rossi X-ray Timing Explorer facility. A possible 
optical counterpart was identified to be a low-mass star (Orosz 
et al. 1998) and the distance to the source was estimated to be 
about 5.3 kpc (Orosz et al. 2002). A variable radio source was 
subsequently found at the optical position (Campbell- Wilson et 
al. 1998) and confirmed by Marshall et al. (1998) based on A5CA 
observations. Radio jets with apparent superluminal velocities 
were observed after the strong X-ray flare in September 1998 
(Hannikainen et al. 2001), placing XTE J 1550-564 among al- 
ready known microquasars. 

The source 2.5-20 keV spectrum was similar to the spectra 
of sources that are dynamically established to be black-holes. 
XTE J 1550-564 has been observed in the very high, high/soft 
and intermediate canonical outburst states of black-hole X-ray 
novae (Sobczaket al. 1999a, 1999b; Cui et al. 1999). The source 
gained additional support for its black-hole nature due to the de- 
tection of a variety of low-/ QPOs (0.08-18 Hz) as well as high- 
/ QPOs (100-285 Hz) during some of the PCA/RXTE point- 
ings (Bradt et al. 1993; Cui et al. 1999; Remillard et al. 1999; 
Wijnands et al. 1999; Homan et al. 2001). 

To probe the object quasi-periodic variability in time-scales 
of seconds we use a pointed observation of the Proportional 
Counter Array detector of RXTE taken from the public archive 
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Fig. 1. (a) An exemplary 100 s data segment extracted from the 
light curve of XTE J 1550-564 as observed by RXTE/PCA on 
September 29th 1998. (b) The histogram of data in our sample 
(1000 s) fitted with a lognormal distribution. 



of HEASARC (|http://heasarc .gsf c .nasa.govll . After Zdziarski & 
Gierlinski (2004), we selected the PCA data set of ObsID of 
30191-01-15-00 containing source observation on 29.09.1998 
(MJD 51085) when the binary system was in its very high state 
(Gierlinski & Done 2003) and displayed narrow QPO feature of 
quahty factor g - 1 1 (Sobczak et al. 2000). 

We reduced the data with the LHEASOFT package ver 6.6.1 
applying the standard data selection: the Earth elevation angle 
> 10°, pointing offset < 0°.01, the time since the peak of the last 
South Atlantic Anomaly SAA > 30 min and the electron con- 
tamination < 0. 1 . The number of Proportional Counting Units 
(PCUs) opened during the observation was equal five. 

We extracted a light curve in 2.03-13.06 keV energy band 
(channels 0-30) with a bin size of Af = 2"^ s for the first part 
of PCA observation (2502 s). For the purposes of this paper, we 
limited our data sample only to the first 1000 s dividing it into 
ten 100 s (A^ = 3200 points) segments for which the wavelet and 
Matching Pursuit analysis was conducted in detail. An exem- 
plary data segment is presented in Fig.[T^. Our sample returned 
the mean countrate and r.m.s. variability equal 1.65 x 10** cts 
s"' and 18.9%, respectively. The data histogram (Fig.[T]3) can be 
fitted with a lognormal distribution. 



f(x;i^,o-) = 



1 



exp 



(- Inx - jj)^ 



2o-2 



(1) 



and the best fit returns the distribution's parameters, the mean 
and standard deviation, equal// = 9.698 + 0.003 and cr = 0.184 + 
0.002, respectively, where the errors are given at 99% confidence 
level (see also Uttley et al. 2005 for a detailed discussion on X- 
ray light curve distribution). 
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Fig. 2. Power spectral density plot calculated based on the entire 
light curve. Main QPO feature is well detected at the frequency 
of 4 Hz. 



We computed Fourier power spectra using the powspec 
subroutine in HEASoft package with Poisson noise level 
subtraction applied. STFT spectrogram plots are derived 
based on the Matlab® Signal Processing Toolbox ver 6.9. 
For the wavelet analysis the Time-Frequency Toolbox 
(|http://tftb.nongnu.org) was used whereas all results corre- 
sponding to the Matching Pursuit signal decomposition were 
computed using MP ver. 4 code provided by Piotr Durka of 
the Department of Biomedical Physics at Warsaw University, 
Poland (http://www.eeg.pl/softwarej. 

3. Time-frequency analysis of 4 Hz QPO structure 

Fig. |2] displays the Fourier power spectral density (PSD) calcu- 
lated from 78 averaged spectra based on 1024 point data seg- 
ments. A QPO peaking at a frequency of 4 Hz is clearly evident 
together with its harmonic structure, and r.m.s. variability calcu- 
lated in a 3.5-5 Hz frequency band equals 12.1%. 

We use the first 100 s of data (Fig.[T^) to compare the three 
different time-frequency techniques. The upper plot in Fig. |3] 
shows a standard STFT spectrogram, computed on a 1 .44 s long 
rectangular (sliding) windo\43. Thus all oscillations of lifetime 
shorter than ~ 1 .44 s are smeared, and the frequency resolution 
is 0.7 Hz, comparable to the FWHM of the QPO. The QPO ap- 
pears to be fairly continuous thoughout the dataset, with small 
excursions around 4 Hz. 

By comparison, the middle panel of Fig.[3]shows the wavelet 
power spectrum using the Morlet wavelet basis functions, i.e. 
modulated sinusoid confined by a Gaussian envelope, (/'(f) = 
;7j--i/4g!(2n-/oOg-' -witii 2nfo - 6 to make each wavelet last for 
about 4 cycles irrespective of the timescale f for the oscillation. 
At the QPO frequency of 4 Hz (wavelet scale a = 0.2427) the 
wavelets can resolve signals which last for 1.44 s. Across the 
QPO width, the wavelet window changes its duration from 1 .63 s 
at 3.5 Hz to 1.14 s at 5 Hz causing a continuous change of fre- 
quency resolution from 0.61 Hz to 0.88 Hz, respectively. This 
gives a better view of the detailed structure of the QPO than the 
STFT spectrogram, but still limits our view of the detailed com- 
position of the QPO signal to half that of the QPO's FWHM. 
Nontheless, this is already able to show that the apparently con- 
tinuous QPO in the STFT result is actually composed of distinct 
events. 



' We select a 1.44 s window size to allow for a better comparison of 
the time-frequency QPO structure resolution between the spectrogram 
and the wavelet analysis at a frequency of 4 Hz. 
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Fig. 3. Time-frequency analysis performed for the XTE J 1550- 
564 light curve presented in Fig.|2^: (upper, a) Spectrogram cal- 
culated for the first 100 s long data segment revealing QPO activ- 
ity around a central frequency of 4 Hz (46 point sliding window); 
(middle, b) the corresponding wavelet power spectrum derived 
with a Morlet analysing function; (bottom, c) Matching Pursuit 
decomposition (based on 3 x 10'' atom dictionary) displayed as 
energy density (Ex){t, f) (Eq. IA.13I) uncovering the main quasi- 
periodic components located mainly between 3.1—^.1 Hz. For aU 
figures colour coding is assumed as follows: increasing values of 
the STFT/wavelet/MP power (energy density) are denoted by a 
gradual brightening of the colour, i.e. from dark blue, yellow to 
red (colour values not in scale). 



Results from the MP algorithm (Section IA.2I ) are shown in 
the bottom panel of Fig. [3] This is deconvolved using 500 it- 
erations and displayed as energy density (Eq. IA.13b . This is 
equivalent to saying that the signal is approximated by 500 sep- 
arate Gabor atoms (hereafter the term atom will be used inter- 
changeably). Together these account for over 95% of total sig- 
nal powefl or speaking conversely, we can reproduce the en- 
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Fig. 4. As in Fig 2, but zooming in on the 3.5-5 Hz frequency 
band dominated by the QPO. 



tire light curve in 95% based on 500 fitted Gabor atoms (figure 
not shown). The energy density has been displayed with a lin- 
ear scale in order to show the time-frequency distribution of the 
strongest events. 

Fig. |4] shows the detailed structure of the QPO (the 3.5- 
5 Hz band) derived from each of these techniques. This imme- 
diately shows the progressive increase in time-frequency reso- 
lution from STFT to wavelet to MP analysis. Since MP algo- 
rithm scans the whole stochastic dictionary for the best match- 
ing Gabor functions to the signal, we gain a new opportunity 
to represent the QPO structure by a number of localised periodic 
functions as defined by Eq. (IA.5b . The MP map shows that much 
of the power of the QPO signal in these data is made from 6 to 
10 periodic oscillations around 4 Hz, with life-times between 
20-60 s. We will return to this result later 

The underlying broadband continuum noise (dominated by 
the intrinsic variability rather than purely Poisson statistics, see 
Fig. 1) is also present in all these plots. This appears in the STFT 



^ Before each MP decomposition we normalise the signal to have unit 
energy. Signal energy carried by a single Gabor atom (also referred to 



as an atom energy) is \(R'x, and the sum over all atoms (Eq. lA.lll 
tends to unity. 
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Fig. 5. The upper panel shows the histogram of atom life-times 
as calculated by the MP decomposition within the 3.5-5 Hz 
band. This clearly peaks at short timescales, but extends con- 
tinuously up to the longest timescales sampled by the data. The 
lower panel shows the power carried by QPO atoms of each 
timescale (see text on atom selection). Unlike the atom lifetimes, 
this does not peak at zero, but instead shows that most of the 
QPO power is carried by oscillations lasting ~5- 10 s. 



as the gradation of background colour from pale to dark from 
the bottom (1 Hz) to the top (10 Hz) of Fig. [3^. In the wavelets, 
this is seen as pale, random smudges around 1-3 Hz (Fig.[3}3), 
whereas in MP they are very short-lived signals (less than 0.5 s) 
forming vertical features across the plot. 



4. Detailed QPO structure from Matching Pursuit 

We can quantify the composition of the QPO by showing the dis- 
tribution of the life-times of 500 Gabor atoms found in the MP 
decomposition accumulated for ten segments of XTE J 1 550-564 
data in 3.5-5 Hz band (Fig.|5^). These show a continuous dis- 
tribution of lifetimes, peaked at zero but extending out to the 
longest timescales sampled by the data. Fig. shows the dis- 
tribution of power carried by QPO atoms of each lifetime (sum 
of atom energy times number of atoms over all the events in the 
time bin of that duration). By QPO atoms we define these events 
distributed between 3.7-4.7 HtQ and exclude these with a dura- 
tion below 0.5 s, or with energy which is less than the events in 
the 2.5-3 Hz bandpass (reference red-noise level; see Fig. [J]) as 
both of these will predominantly trace the broadband noise. 

We found that the largest number of atoms (more than half 
of the total 500 per data segment) are characterised by a time du- 
ration of 0.5-3.5 s, but that these together carry only half the en- 
ergy of the strongest 10 atoms, which have lifetimes of 20-60 s 
(see Fig.Hj;). On the other hand, Fig.|5}' reveals that the major- 
ity of QPO power is carried by QPO events with a duration of 
between 5-10 s. 

We now explore to what level the results can distinguish be- 
tween different models for the broadening of the QPO signal. 



' The frequency interval of a width corresponding to the QPO's full 
width at half maximum. 
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Fig. 6. Simulated QPO signal (a; upper panel) which is fully co- 
herent and single frequency but with random amplitude modu- 
lation. The middle panel (b) shows the wavelet decomposition. 
This is not dissimilar to that from the real data in Fig.|4j3, though 
there are fewer large gaps and the derived frequency is more sta- 
ble. The bottom panel (c) shows the MP decomposition of the 
same signal. This is clearly different to that of the real data in 
Fig. |4j;, ruling out this QPO shape. Colour coding of a wavelet 
and MP map assumed the same as described in a caption of 
Fig.E] 



4. 1. Coherent, single frequency signal with amplitude 
modulation 

Continuous amplitude modulation of a single frequency sinu- 
soid will result in a broad peak in a power spectrum. We sim- 
ulate this by generating a single sinusoid modulation y(t) = 
A sin(27r/of) + C over 100 s with sampling time of At = 2"^ s, 
/o = 4 Hz, and a constant value of C - 1.65 x 10"* cts s"'. 
The amplitude A is taken as a random variable from the log- 
normal distribution as fitted to the XTE J1550-564's histogram 
(see Fig. [TJ5 and Sect. |2]for details). For y(t) we perform its ex- 
ponential transformation, namely exp[37(f)], which accounts for 
a non-linear type of variability present in X-ray light curves of 
accreting black-hole systems (Uttley & McHardy 2001; Uttley, 
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Fig. 7. As in Fig.|6] but for a simulated QPO signal which is a 
fully coherent chirp signal, with frequency increasing with time 
from 3.9 to 4.2 Hz. A smooth frequency drift is readily apparent 
in both wavelet and MP maps, unlike the data shown in Fig.|4j3 
and Hj;. Colour coding of a wavelet and MP map assumed the 
same as described in a caption of Fig. [3] 



McHardy & Vaughan 2005). In this point we ensured the signal 
would have a zero mean and its variance to be such to obtain 
(after the above transformation and Poisson noise inclusion) a 
similar degree of r.m.s. as observed for the entire XTE J1550- 
564 data sample as well as to match r.m.s. variability in the QPO 
band (3.5-5 Hz). 

Fig. |6h displays a simulated signal, whereas Fig. ^ re- 
veals the wavelet power spectrum. The spectrum is similar to 
the XTE J 1550-564 wavelet map, but the MP decomposition is 
clearly different to that in Fig. |4};. The first iteration returned a 
~85 s long Gabor atom exactly at a frequency equal to the as- 
sumed /o, carrying about 85% of the signal energy and dominat- 
ing over the rest of the fitted atoms. 

Therefore the MP decomposition clearly rules out a QPO 
signal formed solely from amplitude modulation as the real data 
cleanly show multiple strong Gabor atoms with frequencies be- 
tween 3.7-4.7 Hz. 



4.2. Coherent chirp signal with amplitude modulation 

Physical models for the QPO involving blobs spiralling inwards 
should produce a monotonically increasing frequency signal, 
otherwise known as a chirp (e.g. Turner et al. 2006; Mhlahlo 
et al. 2007). 

We repeate the experiment above, but with a chirp signal 
where y(t) = A sin[27r/(f)f] + C and / changes continuously be- 
tween 3.9 and 4.2 Hz. Fig. |7] shows that again the amplitude 
modulation can be easily captured in the wavelet map, but it 
leaves no impact on frequency. In MP map, the continuously 
changing frequency is fitted with separate atoms (see Sect. lA.24l 
for more details on Fig. [TJ; construction) which trace the drift. 
The lifetimes of the individual atoms are around 10-20 s and 
show that the data can distinguish between frequencies separated 
by ~ 0.03 - 0.06 Hz, i.e. approximately an order of magnitude 
better than the wavelet technique at these frequencies. 

Again, this MP map is quite unlike that of the real data shown 
in Fig.|4};. 

4.3. Short lifetime, random amplitude signals 

The tests described above clearly show that the QPO signal in 
the data is not given by a single frequency. There are multiple 
frequencies present in the dataset, but these frequencies do not 
show any systematic, long timescale trend such as a chirp. 

This seems to fit most easily into a picture where the QPO is 
composed of separate events with shape rather similar to that of 
the Gabor atom i.e. a single frequency lasting for some timescale 
fqpo which varies. This would fit quite well into models where 
the QPO can be randomly excited by turbulence, and last for 
some timescale before getting damped out. The specific LF QPO 
model of Ingram et al. (2009) has a vertical, Lense-Thirring pre- 
cession of the hot inner flow. This will be damped on a viscous 
timescale, which is of an order ~ 3 s for typical parameters. 

Here, we approximate this physical picture by assuming that 



100 fqpo varies uniformly between (0,3] s and that after each QPO 



event there is a waiting time of fbreak randomly distributed in the 
(0, 2] s interval. Therefore the QPO is qpo(f) = A sin(27r/of + 
(f>) + C where /o - 4Hz and a phase is randomised for every new 
QPO. As previously, we assume the amplitude to be a random 
variable, however this time drawn from the exponential distribu- 
tion. 



100 



/(x;|/) = -exp(--), 



(2) 



fitted to the histogram of amplitudes derived for QPO atoms 
(3.7-4.7 Hz) from the XTE J 1550-564 data, where we obtain 
yii = 2692+^45 at the 99% confidence level (Fig. |9]l. During the 
breaks we approximate signal value by a random variable drawn 
from lognormal distribution and we follow the signal transfor- 
mation as described in Sect. 14.11 

The above construction of simulated signal allows us to ac- 
count for a general countrate distribution to resemble that in 
XTE 1550-564 light curve as well as to approximate more 
adequately the amplitudes of simulated quasi-periodic events. 
Fig. [8^ displays an examplary 100 s signal obtained in the above 
manner During the simulation we assumed the r.m.s. variabil- 
ity in 3.5-5 Hz band to match closely 12% as derived from the 
XTE J1550-564 data. In Fig.[TO]we compare both Fourier power 
spectra, i.e. for the X-ray source and simulated signal (here taken 
to be 1000 s long). The plots indicate on very good qualitative 
agreement between the original and simulated QPO feature in 
the frequency domain, though it is straighforward to note that our 



Lachowicz &. Done: QPOs under wavelet microscope 



7 




20 40 60 80 100 

Time [s] 




20 40 eo SO 100 

Tim 9 [s] 




10 20 3O4O50eO7106O90 10O 
Time [s] 



Fig. 8. As in Fig. |6] but for a simulated QPO signal which is 
composed of separate oscillations which last for a timescale ran- 
domly distributed between 0-3 s followed by a wait time which 
is randomly distributed between 0-2 s. Each oscillation has ran- 
dom phase and amplitude but is at a fixed frequency of 4 Hz. The 
upper panel (a) shows an exemplary light curve. The middle sec- 
tion (b) displays a corresponding wavelet map, while the lowest 
panel (c) shows MP decomposition. Both these are rather similar 
to those from the real data in Fig.|4j3 and|4j;. Colour coding of a 
wavelet and MP map assumed the same as described in a caption 
of Fig. E] 



simulation does not correctly reproduce broad-band red-noise 
variability as observed in XTE J 1550-564. 

Encouraged by the above agreement, we zoom in on the in- 
ternal structure of our simulation by performing wavelet and 
Matching Pursuit analysis as previously (Fig. |8]l, and we find 
that both maps provide the results which look very similar to 
that from the real data (cf. Fig.|33 and c). 

More quantitatively. Fig. [TT| shows that the distributions of 
atom lifetimes and energies are also remarkably well matched 
by the simulation. Even though the signal is made from separate 
oscillations with life-time of less than 3 s, an MP decomposition 
displays a number of longer lasting atoms, e.g. of lifetimes be- 



300 
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Fig. 9. The distribution of amplitudes for QPO Gabor atoms ex- 
tracted from Matching Pursuit analysis of XTE J1550-564 light 
curve in 3.7-4.7 Hz frequency band. The solid line denotes the 
exponential distribution fitted to the data. 
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Fig. 10. The comparison between power spectral density caclu- 
lated for XTE J 1550-564 (open black circles) and simulated sig- 
nal as described in Sect. l4.3l (filled red circles). 

tween 40-80 s. Since this issue requires a better understanding, 
we examined the MP response to the set of various test signals 
containing different noise levels, amplitude distribution, variable 
frequency and/or constant phase for separate qpo(f) events. In 
conclusion, we have noticed that the longer lasting atoms can 
be fitted to the data by an MP algorithm making use of an atom 
concatenation process. In this case, a number of shorter lasting 
oscillations characterized by very similar phase and frequency 
(not amplitude!) can be linked together into one longer atom. 
Therefore, e.g. a 60 s long atom can describe, say, three to five 
~3 s oscillations, and our inability to distinguish among them 
refers to the limited time-frequency resolution of the Matching 
Pursuit algorithm. 

Similarly, while the simulation has only one constant fre- 
quency signal, there is a spread in frequency in the MP decom- 
position which matches that seen in the data. 

The low-frequency QPO seen in the real data can be well 
described by a signal which consists of random, discrete events 
at a single frequency, with a spread in lifetime and amplitude. 

5. Discussion and Conclusion 

Quasi-periodic oscillations observed in the power spectra of 
many accreting compact objects still remain the most intriguing 
puzzle. Since high-frequency QPOs are generally thought as the 
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Fig. 11. As in Fig. |5] histograms of atom lifetimes and energies 
for the simulated signal as described in Sect. 14.31 Bottom panel 
refers to QPO atoms defined in 3.7-4.7 Hz band, excluding these 
events of duration less than 0.5 s and of atom energy lower than 
in 2-3.7 Hz and 4.7-5 Hz bands. 



most attractive landmark feature used in testing the behaviour 
of matter in the very strong gravitational field (e.g. Kluzniak & 
Abramowicz 2001; Rebusco & Abramowicz 2006), on the con- 
trary, the low-frequency QPOs have been found to be more enig- 
matic in their nature (Casella et al. 2006). For the latter, a quanti- 
tative analysis has been a dominant tool towards our comprehen- 
sion of the global QPO properties (e.g. Cui et al. 1999; Sobczak 
et al. 2000). 

The presence of broad peaks in the Fourier power spectra fit- 
ted with Lorentzian profiles suggests a possible mechanism of 
X-ray engine at work: damping of locally emitted radiation at 
the very short time-scales (e.g. Misra & Zdziarski 2008). The 
position of QPO peak frequency varies in time, generally being 
drifting between 0.1-10 Hz (van der KUs 2005, 2006). It has 
been suggested it may be correlated with both underlying geom- 
etry of the inner flow, and accretion rate, thus with the spectral 
state of the source (e.g. Di Matteo & Psaltis 1999; Dutta et al. 
2008; Qu et al. 2010). The r.m.s. variabiUty of LF QPOs is also 
significant and in some cases exceeds 10%. 

The successful modelling of this QPO phenomena requires a 
sophisticated approach which would include acceptable model's 
results not only in the reproductions of the light curve and 
Fourier PSD, but also in providing an agreement in a func- 
tion of energy (e.g. Sobolewska Zycki 2006). Such attempts 
have been so far conducted on many occassions (e.g. Giannios 
& Spruit 2004; Machida et al. 2006; Machida & Matsumoto 
2008). Fundamental and dynamical Fourier analysis of X-ray 
light curves gained a leading position in spoting the existence 
and evolution of QPOs (e.g. Galleani et al. 2001; Pottschmidt et 
al. 2003; Barret et al. 2005) but failed in delivering the satisfac- 
tory answer on the physical structure due to insuflicient resolu- 
tion. On the hand, the application of new time-frequency meth- 
ods for X-ray time-series analysis (e.g. by wavelet approach) 
initiates a novel look at the data and brings an opportunity for 
testing QPO models in a new space of parameters (e.g. Steiman- 
Cameron et al. 1997; Pai'ker 1998; Lachowicz & Czerny 2005). 

In this paper we have applied the Matching Pursuit method 
for signal decomposition and shown how MP technique can give 
much better resolution in the time-frequency plane to study the 
evolution of the low-frequency QPO seen from black hole binary 
systems. We have used data from the bright transient systems 



XTE J 1550-564 to illustrate this with an observation taken in 
the very high state (aka the steep power law state) where the 
QPO is strong and relatively coherent. 

The Mathching Pursuit decomposition clearly rules out sev- 
eral potential origins for the broadening of the QPO signal. 
The 'quasi' nature cannot be due to a single frequency, coher- 
ent oscillation with amplitude modulation, nor can it be a long 
timescale (~100 s) systematic change in frequency (chirp sig- 
nal). Instead, it is well matched by a series of discrete oscilla- 
tions, with life-times between 0-3 s, each with random phase 
and amplitude but fixed frequency. 

Physically, this is consistent with models of the QPO where 
it is excited by random turbulence. In particular, the Lense- 
Thirring precession model of Ingram et al. (2009) can fit nicely 
into this picture, where the vertical precession is excited by the 
large scale stochastic turbulence generated in the accretion flow. 
The timescale for decay of the vertical precession should be re- 
lated to the viscous timescale of the flow. In the Ingram et al. 
(2009) model, a 4 Hz QPO would be associated with an outer 
radius of the hot flow of ~ lORg, so should have a viscous 
timescale of ~ 2 - 3 s for the typical parameters used by Ingram 
et al. (2009). 

The application of Matching Pursuit opens a brand new op- 
portunity for studies of LF QPOs in the time-frequency domain. 
It encourages for revision of archival data of other BHB which 
display narrow QPOs in their Fourier power spectra. By in- 
creased time-frequency resolution provided by MP it is possible 
to diferentiate amongst various QPO models, as demonstrated in 
this paper Therefore, the Matching Pursuit analysis may bring 
us much closer towards understanding of the underlying physical 
properties of the QPO phenomenon. 
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Appendix A: Time-Frequency Analysis 

A.1. Spectrogram and wavelet analysis 

In X-ray time analysis the wavelet transformation constitutes an 
alternative tool for signal decomposition when confronted with 
the classical approach based on the Fourier transform (PSD). In 
the latter, a signal is analysed with the basis functions ^ oc e'-''. It 
detects strictly periodic modulations (with period of 1 //) buried 
in the signal. Therefore, no information about time-frequency 
properties is available. 

A simple solution to this problem has been found via the ap- 
plication of Short-Term Fourier Transform (STFT). Basic con- 
cept standing for this method is to divide signal into small seg- 
ments and calculate for each of them Fourier transform. The 
signal analysis performed using spectrogram, i.e. |STFTp often 
blurs all signal oscillations of life-times less than the duration of 
a segment causing a loss of significant information. 

The wavelet analysis treats this problem more accurately as a 
basic concept standing for its usefulness relay on probing time- 
frequency plane at different frequencies with different resolu- 
tions. The wavelet power spectrum for a discrete signal can be 
defined as the normalised square of the modulus of the wavelet 
transform ^|w„(a,„)p where ^ stands for normalisation factor 
and w„(fl,„) a discrete form of the continuous wavelet transform 
(Farge 1992; Torrence & Compo 1998; Lachowicz & Czerny 
2005). It is attractive to adopt a Morlet wavelet to probe quasi- 
periodic modulations. Since Morlet wavelet oscillates due to a 
term oc e", it is perfectly suited for this task. 

A.2. Matching Pursuit algorithm 
A.2.1. Parlez-vous frangais? 

Signal time-frequency analysis can be compared to speaking in 
a foreign language. In each language we use words. Words are 
needed to express our thoughts, problems, ideas, etc. By a smart 
selection of proper words we can say and explain whatever we 
wish. A whole collection of words can be gathered in the form 
of a dictionary. One can express simple thoughts using a very 
limited set of words from a huge dictionary (a subset). The same 
can be applied to a time-series analysis. In order to describe the 
signal one needs to use a minimum available set of functions - 
orthonormal basis functions. 

Description of the signal x„ which uses only basis functions 
g„ in - \, ...,N') is therefore limited. This is often the case with 
wavelet transform. In order to improve signal description, one 
can increase the size of a basic dictionary by adding extra func- 
tions. Such a dictionary's enlargement introduces a redundancy 
what is a natural property of all, for instance, foreign languages. 
The most reliable signal description can be achieved by a signal 
approximation: 

M<N' 

X * ^ ay^ngy.n (A.l) 
n=l 

where coefficients defined simply as the inner products 

of basis functions ^y_„ with a signal: 

Xco M<N' 
X{t)gy{t)dt X ^ X„gy^„ (A.2) 
n=\ 

and y refers to a set of parameters of function g. 

In each approximation, a relatively small number of func- 
tions is generally required. For each case the signal is projected 
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onto a set of basis functions, therefore, in other words, its rep- 
resentation is always fixed. However, what if one may approach 
the signal decomposition problem from the other side, i.e. to fit 
the singal representation to the signal itself? One can achieve it 
by choosing from a huge dictionary of pre-defined functions a 
subset of functions which would best represent the signal. Such 
solution had been proposed by Mallat & Zhang (1993) and is 
known under the name of Matching Pursuit algorithm. 

A.2.2. Algorithm 

Matching Pursuit (MP) is an iterative procedure which allows 
to choose from a given, redundant dictionary D, where D - 
{gu82,-,8n] such that \\gi\\ = 1, a set of m functions {g;),„<„ 
which match the signal as well as possible. D is defined as a 
family (not a basis) of time-frequency waveforms that can be 
obtained by time-shifting, scaling and, what is new compar- 
ing to wavelet transform, modulating of a single even function 
geLHn)-\\g\\^l- 



1 

ya 



t-b 



Jf 



where an index y refers to a set of parameters: 



(A.3) 



(A.4) 



Each element of the dictionary is called an atom and usually is 
defined by a Gabor function, i.e. a Gaussian function modulated 
by a cosine of frequency /: 



gy(t)^Ky(t)Gy(t)^Ky(t)eXp 





it-b\'] 








cos 




[ a j\ 





COS (27r/f -H 0) (A.5) 



where 



KM) = 



2 



and a denotes the atom's scale (duration; in this paper also re- 
ferred to as atom's life-time) and b - time position. A selection 
of Gabor functions for purposes of MP method is fairly done 
as they provide optimal joint time-frequency localisation (e.g. 
Mallat & Zhang 1993; Cohen 1995; Flandrin 1999). 

In the first step of MP algorithm, a function gy i is chosen 
from D which matches the signal x best. In practice it simply 
means that a whole dictionary D is scanned for such gy^i for 
which a scalar product \{x,gy^i}\ e is as large as possible. 
Therefore, a signal can be decomposed into: 



X = {X,gy,l}gy,l +R^X 



(A.7) 



where a a residual signal, has been denoted. In the next step, 
instead of x, the remaining signal is decomposed by finding 
a new function, gy^, which matches R^x best. Such consecu- 
tive steps can be repeated p times where /? is a given, maximum 
number of iterations for signal decomposition. In short, MP pro- 
cedure can be denoted as: 



R^h = X 

R'x^{R'x,gyj)gyj + R'^'x 
gyj = argmaXj^^,eo \{R'x,gyj,)\. 



(A.8) 



Davis et al. (1994) showed that for x e H the residual term R'x 
defined by the induction equation ( IA.8b satisfies: 



lim \\R'x\\ = 



(A.9) 



i.e., when H is of finite dimension, ||^'x|| decays to zero. Hence, 
for a complete dictionary D, the Matching Pursuit procedure 
converges to x: 



Y^{R'x,gyj}gyJ 



(A. 10) 



1=0 



where orthogonality of R'^ x and gyj in each step impUes energy 
conservation: 



(A. 11) 



A.2.3. Time-frequency representation 

From (lA.lOb one can derive a time-frequency distribution of sig- 
nal's energy by adding a Wigner-Ville distribution of selected 
atoms gy: 

(Wx)(t, f) ^ZZo 8y,d\Hwgy,d(t, /)+ 

ZSo 27=0 J^i(R'^^ gyjKR'^, gy.,T X (A. 12) 

X mgyMgyjWJ) 

The double sum present in the above equation corresponds to 
cross-terms present in the Wigner-Ville distribution. In MP ap- 
proach, one allows to remove these terms directly in order to 
obtain a clear picture of signal energy distribution in the time- 
frequency plane (f, /). Thanks to that, energy density of x can be 
defined as follows: 



{Ex){t,f) = £ \{R'x,gy,i)\\WgyMUD- 



(A. 13) 



Since Wigner-Ville distribution of a single time-frequency atom 
gy satisfies: 



(Wgy,d{tJ)dtdf^\\gy\\^^\ 

OO %J —oo 



(A. 14) 



thus combining < IA.14I > with energy conservation of the MP ex- 
pansion < IA.lQl t yields: 



oo %J —i 



iEx)it,f)dtdf^\\x\\ 



(A. 15) 



what justifies interpretation of (Ex)(t, f) as the energy density of 
signal X in the time-frequency plane. 

In practical computation of MP algorithm, sampling of a dis- 
crete signal of length N - 2^ is governed by an additional pa- 
rameter j: an octave (Mallat & Zhang 1993). The scale a corre- 
sponding to an atom's scale is derived from a dyadic sequence 
of a = 2 ', < j < L. Thus, parameters b and /, corresponding to 
atom's position in time and frequency, respectively. By resolu- 
tion in the MP one may understand the distance between centres 
of atoms neighboring in time or in frequency. A dyadic sam- 
pling pattern proposed by Mallat & Zhang has been found to suf- 
fer from a statistical bias introduced in the atom parametrisation 
(IA.4b . Durka et al. (2001) proposed a solution to this problem by 
application of stochastic dictionaries. In this approach the pa- 
rameters of dictionary's waveforms are randomised before each 
decomposition by drawing their values from continuous ranges. 

In this paper in all time-series analyses with MP method, we 
decompose signal in 500 iterations using, if not otherwise stated, 
stochastic dictionaries composed of 3 x 10^ atoms. By the term 
strongest atoms we will refer to these Gabor atoms for which 
the calculated inner product with a signal has been largest (in 
practice it refers to the atoms fitted within first 1-10 iterations). 
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A.2.4. Number of MP decompositions for a single signal 

Durka et al. (2001) pointed at a significant aspect of Matching 
Pursuit method: an ability to represent signals with continuously 
changing frequency. In the standard approach one aims at cal- 
culation of a single MP decomposition for a given signal us- 
ing a large redundant dictionary of pre-defined functions (e.g. 
3 X 10^ atoms). When the signal contains continuously changing 
frequency in time a single MP decomposition will provide en- 
ergy density plot to be composed of a number of single Gabor 
atoms. Therefore, such a view remains far different from expec- 
tations (see exemplary illustration in Durka et al. 2001, 2007). 

On the way of experiments Durka et al. (2001) discovered 
that by introduction of stochastic dictionaries and averaging a 
larger number, say 100, of the MP energy density plots together 
one can uncover continuously changing signal frequency. The 
trick is to use for every MP decomposition performed for the 
same signal a new stochastic dictionary and of much smaller 
size, e.g. composed of 5 x lO** atoms only. It has been applied by 
us in order to create Fig.|7};. 

In this paper we have checked that for all our data segments 
the MP results are insensitive to the way of MP computation, 
i.e. irrespective of the dictionary size used and the application 
of aforementioned MP map averaging process. In this way, we 
are able to exclude a scenario in which QPOs are being gener- 
ated due to long timescale oscillations of continuously changing 
frequency (see Section l4r2l for details). 



